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We simulate late-stage coarsening of a 3D symmetric binary fluid using a lattice Boltzmann method. 
With reduced lengths and times, I and t respectively (with scales set by viscosity, density and 
surface tension) our data sets cover 1 H I H 10^, 10 ^ t J; 10*. We achieve Reynolds numbers 
approaching 350. At Re ~ 100 we find clear evidence of Furukawa's inertial scaling {I ~ t^ ), 
although the crossover from the viscous regime [l ~ t) is very broad. Though it cannot be ruled out, 
we find no indication that Re is self-limiting {I ~ i^") as proposed by M. Grant and K. R. Elder 
[Phys. Rev. Lett. 82, 14 (1999)]. 
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When an incompressible binary fluid mixture is 
quenched far below its spinodal temperature, it will 
phase separate into domains of different composition. 
Here we consider only fully symmetric 50/50 mixtures 
in three dimensions, for which these domains will, at 
late times, form a bicontinuous structure, with sharp, 
well-developed interfaces. The late-time evolution of this 
structure remains incompletely understood des pite theo- 
retical Q-Ql , experimental |^ and simulation p- flQ] work. 

As emphasized by Siggia ||l| and Furukawa^, the 
physics of spinodal decomposition involves capillary 
forces, viscous dissipation, and fluid inertia. Thus, as- 
suming that no other physics enters, the control param- 
eters are interfacial tension a, fluid mass density p, and 
shear viscosity 77. From these can be constructed only 
one length, Lq = 'ff' / P^^ and one time Tq = rf / pa^. 
We define the lengthscale L{T) of the domain struc- 
ture at time T via the structure factor S{k) as L — 
2tt j S{k)dk/ J kS{k)dk. The exclusion of other physics 
in late stage growth then leads us to the dynamical scal- 
ing hypothesis ||,|^: I — l{t), where we use reduced time 
and length variables, I = i/Lo and t = {T — Tint)/TQ. 
Since dynamical scaling should hold only after interfaces 
have become sharp, and transport by molecular diffu- 
sion ignorable, we have allowed for a nonuniversal offset 
Tint', thereafter the scaling function l{t) should approach 
a universal form, the same for all (fully symmetric, deep- 
quenched, incompressible) binary fluid mixtures. 

It was argued further by Furukawa 0] that, for small 
enough t, fluid inertia is negligible compared to viscosity, 
whereas for large enough t the reverse is true. Dimen- 
sional analysis then requires the following asymptotes: 

l^bt; t-^t* (1) 

t » t* (2) 
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where, if dynamical scaling holds, amplitudes 6, c (and 
the crossover time t*) are universal. The Reynolds num- 
ber, conventionally defined as, Re = p/rjLdL/dT — II, 
becomes indefinitely large in the inertial regime, Eq. (0). 



In a recent paper. Grant and Elder have argued H] 
that the Reynolds number cannot, in fact, continue to 
grow indefinitely. If so, Eq. (0) is not truly the large t 
asymptote, which must instead have I ^ i" with a < 2- 
Grant and Elder argue that at large enough Re, turbu- 
lent remixing of the interface will limit the coarsening 
rate Q, so that Re stays bounded. A saturating Re 
(which they estimate as Re ~ 10 — 100) would require 
any t^/^ regime to eventually cross over to a limiting t^/^ 
law. But if a single length scale I ~ i^^^ is involved, a 
saturating Re implies balance between viscous and in- 
ertial terms (t~'^"), while the driving term (interfacial 
tension) remains much larger than either (t^^). This 
suggests a failure of scaling altogether, with at least two 
length scales relevant at late times. In any case, the ar- 
guments of Grant and Elder are far from rigorous; the 
coarsening interfaces could, remain one step ahead of the 
remixing despite an ever-increasing Re which, if applied 
to a static interfacial structure, would break it up. Thus 
Eq. (y) cannot yet be ruled out as a limiting law. 

In what follows we present the first large-scale simu- 
lations of 3D spinodal decomposition to unambiguously 
attain a regime in which inertial forces dominate over 
viscous ones. We find direct evidence for Furukawa's 
I ~ t^^^ scaling, Eq. (||). Although a further crossover to 
a regime of saturating Re cannot be ruled out, we find 
no evidence for this up to Re ~ 350. Our work, which is 
of unprecedented scope, also probes the viscous scaling 
regime [Eq. dy)] , and the nature of the crossover between 
this and Eq. (0) . Full details of our results |ll[ and of the 
simulation algorithm p2l will be published elsewhere. 

Our simulations use a lattice Boltzmann (LB) method 
||l3| , |l4| with the following model free energy: 

F = jdv |-|02 + |04 + ^1^^^ ^i^^pj (3) 

in which A, B and k are parameters that determine 
quench-depth {A/ B — > 1 for a deep quench) and inter- 
facial tension (u = ^J%kA^ /^B'^); (/) is the usual order 



parameter (the normalized difference in number density 
of the two fluid species); p is the total fluid density, which 
remains (virtually) constant throughout fl4||l5| |. 

The simulation code follows closely that of ||lj| (for 
details see |0,n2|) and uses a cubic lattice with nearest 
and next-nearest neighbor interactions (D3Q15). ft was 
run on Cray T3D and flitachi SR-2201 parallel machines 
with system sizes up to 256'^. The LB method allows 
the user to choose r/, a, p (we set p = f without loss of 
generality), along with the order-parameter mobility M 
defined via </) = V-M'V{6F/S(j)). Although it plays no role 
in the arguments leading to Eqs. (|y) and (H), M must be 
chosen with some care to ensure that at late times (a) the 
algorithm remains stable, (b) the local interfacial profiles 
remain close to equilibrium (so that a is well-defined), 
and (c) the direct contribution of diffusion to coarsening 
is negligibly small. Table | shows the parameters used 
for our eight 256'^ runs. 

TABLE L Parameters used in LB runs 



Lo 


To 


A,B 


K 


V 


M 


O" 


36 


930 


0.083 


0.053 


1.41 


0.1 


0.055 


5.9 


71 


0.0625 


0.04 


0.5 


0.5 


0.042 


5.9 


71 


0.0625 


0.04 


0.5 


0.2 


0.042 


0.95 


4.5 


0.0625 


0.04 


0.2 


0.3 


0.042 


0.15 


0.89 


0.00625 


0.004 


0.025 


4.0 


0.0042 


0.010 


0.016 


0.00625 


0.004 


0.0065 


2.5 


0.0042 


0.00095 


0.00064 


0.00313 


0.002 


0.0014 


8.0 


0.0021 


0.00030 


0.00019 


0.00125 


0.0008 


0.0005 


10.0 


0.00082 



In all runs, the interface width is ^ ~ 5y'K/2A ~ 3 in 
lattice units. This was found |ll| to be the minimum ac- 
ceptable to obtain an accurately isotropic surface tension. 
To minimise diffusive effects, data for which the diffusive 
contribution to the growth rate was greater than 2% was 
discarded |16|| ; this corresponded to a minimum value of 
L of 15 < Lmin < 24, depending on the run parameters. 
The large size of our runs allowed a ruthless attitude to 
finite size effects: we use no data with L > A/4, with 
A the linear system size. In our 256'^ runs, these filters 
mean that the good data from any single run lies within 
20 ~ L < 64, a comparable range to previous studies 
[^-||. Datasets of high and low Lq are well fit respec- 
tively by a = 1 and a — 2/3 (see Fig.|]). 

However, as emphasized by Jury et al. |§], meaningful 
tests of scaling are best made not by looking at single 
data sets but by combining those of different parameter 
values. To this end, the good data from each run were 
fit to L = i3(T — Tint)", so as to extract an intercept 
Tint', we then transformed the data to reduced physical 
units I and t defined above. The exponent a was first 
allowed to float freely; this gave reproducible values at 
large I and t (e.g., a — 0.69 and 0.67 for the last two 
data sets in Table 1), but more scattered ones at small 
I and t {a = 0.88, 0.86 and 1.16 for the first three data 




8000 10000 12000 14000 16000 18000 
T (time steps) 

FIG. 1. L vs. T for the runs shown in Table Q with Lq — 
5.9 (M = 0.2) (circles) and 0.0003 (diamonds). The region 
used for fitting is delimited by {Lmin < L < Lmax ~ 64) and 
the fits are projected back to show the intercepts, Tint- The 
fits (solid) are to a = 1, 2/3; free exponent fits are also shown 
(dashed), with best fit values a = 1.16, 0.69. 



sets). In the latter region the floating fit is relatively 
poorly conditioned; it also gives large relative errors in 
Tint (see Fig.l). In contrast, fits to a = 1 for these three 
data sets gave much better data collapse with consistent 
values oib {b = 0.073, 0.073 and 0.072 ± 0.01). Thus we 
are confident of a = 1 in this region. For the remaining 
data sets we estimate errors in individual exponent values 
at around 10% and in reduced time t around 3% to 10%. 
Figure || shows all our data sets on a single plot using 
reduced variables I and t. Such a plot is necessarily log- 
log, since our data sets span seven decades in t and five 
in ^, a range which exceeds all previous studies combined. 

These LB results are fully consistent with the existence 
of a single underlying scaling curve I — l{t), in which 
viscous {I = bt) and inertial (I = ct^/^) asymptotes are 
connected by a long crossover whose breadth justifies our 
use of a single floating exponent a in the fits used above 
to extract Tint for each run. Although we cannot rule 
it out for still larger times t, we see no evidence for a 
further crossover to a regime with asymptotic exponent 
a < 1/2 as demanded by Grant and Elder |Q. 

Before considering our results in more detail, we dis- 
cuss their relation to others previously published. We 
restrict attention to those 3D data sets for which reli- 
able estimates of Lq and Tq exist ||]. Datasets of Laradji 
et al. M and of Bastea and Lebowitz H are shown on 
Fig.0 (fitted to a = 1 g). These lie in an l,t range 
(1 $/ ^ 20) in which our own data shows viscous (linear) 
scaling [Eq. (||)]; both data sets were claimed to confirm 
the linear law by their authors, but with differing values 
of 6 = 0.13, 0.3. Our own b values are lower than either 
(see above and Fig.||). As noted above, we took special 
care to ensure that the diffusive contribution to coarsen- 
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FIG. 2. Scaling plot in reduced variables (I, t) for LB data, 
bold lines (left to right) are those of Table I (top to bottom). 
Squares Ref. [ho|, triangles Ref. |6|, circles Ref. JtI. Inset: 
DPD data of Ref. |8| (solid lines) with one of our data sets 
(I/O = 0.15, circles) repeated for comparison. 

ing was small; we have found that, for matching Lq,To 
values, LB data sets similar to those of Refs. |gjj can be 
generated using too large a mobiUty M. We hypothesize 
therefore that both data sets have strong residual diffu- 
sion, leading to an overestimate of b. Likewise the data 
of Appert et al. |10) , which lies in the crossover regime 
of our scaling plot, asymptotes to our data from above; 
this suggests that their fitted exponent a ~ 2/3 is too 
low because of diffusion. 

A different explanation, based on a possible nonuni- 
versality of the physics of topological reconncction of do- 
mains, was suggested by Jury et al. ^, whose dissipative 
particle dynamics (DPD) results also appear in Fig.0 (in- 
set) |1^. These authors found that each data set was 
well fit by a linear scaling, Eq. (|l|), but with a system- 
atic increase of the b coefficient upon moving from upper 
right to lower left in the scaling plot [gj. Their alter- 
native suggestion was that their own data, and that of 
Refs. @,0l, were part of an extremely broad crossover re- 
gion, 1 ^ t 1i 10** in reduced time. Our LB data support 
the idea of a broad crossover, but instead places it at 
10^ < i < 10^. Note that, unhke those of Refs. [|,0, 
all the data sets of Jury et al. do lie very close to our 
own (Fig. inset). Since the two simulation methods are 
entirely different, this lends support to the idea of a uni- 
versal scaling, although the fact that each DPD run is 
best fit by a locally linear growth law does not B. The 
latter could be partly due to finite size effects; to obtain 
enough data. Jury et al. included results up to L = A/2, 
whereas we reject all data with L > A/4. 

The arguments of Ref. M involve the intrusion of a sec- 
ond length scale, alongside Lq, which in the LB context 
is the interfacial width ^ (or more generally, a molecular 
scale). The ratio h = ^/Lq for real fluids is in the range 



0.05 (water) to 10~^ (glycerol). In simulations, ^ cannot 
be smaller than the lattice spacing, and the inertial region 
is achieved by setting Lq <^ 1, so h :^ 1. In this sense 
our interface is "unnaturally thick" : simulation runs that 
enter the inertial regime do so directly from a diffusive 
one, without an intervening viscous regime. However, 
this should not matter if l(t) follows a universal curve, 
as our results (in contrast to Ref. ||]), in fact suggest. 
But the microscopic length still plays an interesting role, 
as follows. As a fluid neck stretches thinner and thinner 
before breaking, it shrinks laterally to the scale ^; diffu- 
sion then takes over to finish the job of reconncction. So, 
although our work involves length scales where the direct 
contribution of diffusion to domain growth is negligible, 
we must ensure that it is handled correctly at smaller 
scales. This factor limits the accessible range of / and t, 
not only at the lower [g| but also at the upper end |ll[ . 

The breadth of the viscous-inertial crossover is some- 
what less extreme when expressed in terms of Re (see 
above); our data span 0.1 S: Re ^ 350 and the crossover 
region is roughly 1 S Re S 50. Re values (at L ~ 50) 
for each run are shown in FigJS against reduced time t. 
Data are consistent with Re ~ V^^ as predicted from Eq. 
(0). Note that, in simulating high Re flows, one should 
strive to ensure that the dissipation scale Q (defined as 
^d = {"n^ I^P^Y J with e the energy dissipation per unit 
volume) always remains larger than the lattice spacing. 
This ensures that any turbulent cascade (whose shortest 
scale is A^) remains fully resolved by the grid. Equat- 
ing dissipation with the loss of interfacial energy, one has 
e ~ d{(j/L)/dT and so, in reduced units, A^ — (I? jlf-^^- 
Comparable e values are found directly from our simu- 
lated velocity data; and A^ remains larger than the grid 
size for all our runs p9[ . 

A decisive check that we really are simulating a regime 
where inertial forces dominate over viscous ones, is based 
directly on the velocity fields found in our simulations 
|l]| . From these we calculated rms values of the in- 
dividual terms in the Navier-Stokes equation (p — 1), 
{dv/dt -\- V • Vv) = rjSI^-v — V ■ P. Here P, the pressure 
tensor, contains the driving terms arising from interfa- 
cial tension. Ratios i?i = {dv/dt)rms/{v'^^'^)rms and 
i?2 = (v-Vv)rms/(?7V^v)rms, wcre then computed; these 
can be seen in Fig||. The ratio i?2 is closely related to 
the Reynolds number Re: it differs in representing length 
and velocity measures based on the rms fiuid ffow rather 
than on the interface dynamics and, because the length 
scales associated with the velocity gradients are smaller 
than the domain size, is significantly smaller than Re. 
The dominance (by a factor ten) of inertial over viscous 
forces is, at late times, nonetheless clear (Fig.g). 

We finally ask whether, at the largest Re values we 
can reach, there is in fact significant turbulence in the 
fluid flow. One quantitative signature of turbulence is 
the skewness S of the longitudinal velocity derivatives; 
this is close to zero in laminar flow but approaches S — 
—0.5 in fully developed turbulence |jl^. We do detect 
increasingly negative S as Re is increased but reach only 
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FIG. 3. Reynolds numbers Re — II 
L = 50 for (left to right) runs in Table 
Ratios _Ri (circles), R2 (triangles), the rms inertial to viscous 
ratios (see text) at L = 30 for runs with (left to right) Lq = 
36, 2.9, 0.59, 0.15, 0.054, 0.024, 0.01, 0.01, 0.0016, 0.00095, 
0.00039, 0.0003 (system sizes 96^ (open) and 128^ (filled)). 
Errors are of the order of the symbol size. 

S ~ -0.3 for Re ~ 350 fnt]. This suggests that at our 
highest Re's, turbulence is at most partially developed - a 
view confirmed by visual inspection of velocity maps pi[ . 
Grant and Elder's suggestion of an eventual transition to 
turbulent remixing thus remains open. 

In conclusion, we have presented LB simulation data 
for 3D spinodal decomposition which spans an unprece- 
dented range of reduced time and length scales. At 
t ^ 10^ (Re ~ 1) we observe linear scaling, as announced 
in the previous literature ^^ . This is followed by a long 
crossover (10^ i: t ^ 10^, or 1 i; Re ^ 50) connecting to 
a regime in which inertial forces clearly dominate over 
viscous ones (see FigJs); our work is the first to unam- 
biguously probe this regime in 3D [ |lO| . In the region so 
far accessible (10*^ < t < 10^ or 50 < Re < 350) Furukawa's 
prediction of t'^^^ scaling is obeyed, to within simulation 
error. An open issue is whether this regime marks the 
final asymptote or whether a further crossover occurs to 
a turbulent remixing regime (saturating Re) as proposed 
by Grant and Elder B. If it does, we have shown that 
any limiting value of Re must significantly exceed their 
estimate of 10 — 100. 

We thank Craig Johnston, Simon Jury, David Mc- 
Comb, Patrick Warren and Julia Yeomans for valuable 
discussions. Work funded in part under the EPSRC E7 
Grand Challenge. 
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